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^*^ In recent years there have been many proposals as flexible alternatives to Gaussian based contin- 

uous time stochastic volatility models. A great deal of these models employ positive Levy processes. 
tn , Among these are the attractive non-Gaussian positive Ornstein-Uhlenbeck (OU) processes proposed by 

^^ ' Barndorff-Nielsen and Shephard (BNS) in a series of papers. One current problem of these approaches 

is the unavailability of a tractable likelihood based statistical analysis for the returns of financial as- 
sets. This paper, while focusing on the BNS models, develops general theory for the implementation of 
statistical inference for a host of models. Specifically we show how to reduce the infinite-dimensional 
process based models to finite, albeit high, dimensional ones. Inference can then be based on Monte 
Carlo methods. As highlights, specific to BNS we show that an OU process driven by an infinite activity 
^Nl ' Gamma process, that is an OU-F, exhibits unique features which allows one to exactly sample from 

^ ^ relevant joint distributions. We show that this is a consequence of the OU structure and the unique 

calculus of Gamma and Dirichlet processes. Owing to another connection between Gamma/Dirichlet 
processes and the theory of Generalized Gamma Convolutions (GGC) we identify a largo class of mod- 
_j. ' els, we call (FGGC), where one can perfectly sample marginal distributions relevant to option pricing 

f*^ ^ and Monte Carlo likelihood analysis. This involves a curious result, we establish as Theorem 6.1. We also 

>D ' discuss analytic techniques and candidate densities for Monte-Carlo procedures which can be applied to 

^i^* more general classes of models. 
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1 Introduction 

Barndorff-Nielsen and Shephard (2001a, b)(BNS) introduce a class of continuous time stochastic 
volatility (SV) models that allows for more flexibility over Gaussian based models such as the Black- 
Scholes model[see Black and Scholes (1973) and Merton (1973)]. Their proposed SV model is based 
on the following differential equation, 

(1) dx*{t) = {n + I3v{t))dt + v^l"^ {t)dw{t) 

where x*{t) denotes the log-price level, w{t) is Brownian motion, and independent of w{t), v{t) 
is a stationary Non-Gaussian Ornstein-Uhlenbeck (OU) process which models the instantaneous 
volatility. This latter point is equivalent to the fact that for A > 0, 

v{t)=c-^'v{{))+c~^' [ e^yZ{d\y) 
Jo 

and arises as the solution of the following differential equation, 

dv{t) = -Xv{t) + dZ(Xt). 
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In the above framework Z is a positive homogeneous process, otherwise known as a subordinator, on 
[0,00) and v{0) is an arbitrary positive random variable independent of Z. That is Z(t) := L Z{dy) 
is a stationary process, with Z{Q) = 0, and its distribution specified by its Laplace transform for 
each a; > 0, 

(2) E[e-"^^(*)] = e-*'^(") 

where '0(^) — In (1 ~ c^^'^)p((is), is often called the Levy exponent of an infinite divisible random 
variable equivalent in distribution to Z{1), and p is its corresponding Levy density. Either of these 
characterizes the distribution of the process Z. Importantly, it is obvious from (2), that one does 
not need explicit knowledge of p to calculate ip. Note further that if we wish v{t) to be stationary 

it is necessary to choose v{0) = J_ c^Z*{ds), where Z* is independent of Z but otherwise has the 
same law. 

The model described above is an extension of the Black-Scholes or Samuelson model which arises 
by replacing v with a fixed variance, say cr^. The additional innovation in BNS is that modeling 
volatility as a random process, v{t), rather than a random variable, not only allows for heavy- 
tailed models, but additionally induces serial dependence. This serial dependence is used to account 
for a clustering affect referred to as volatility persistence. The work of Carr, Geman, Madan, and 
Yor (2003) discuss this point further. See also Duan (1995) and Engle (1982) for different approaches 
to this type of phenomenon. The model of BNS has gained a great deal of interest with some related 
works including Carr, Geman, Madan, and Yor (2003), Barndorff- Nielsen and Shephard (2003), 
Eberlein (2001), Nicolato and Venardos (2001), Benth, Karlsen, and Reikvam (2003). See also the 
discussion section in Barndorff-Nielsen and Shephard (2001a). See Carr and Wu (2004) and Duffie, 
Pan and Singleton (2000) for many other models. 

Note that the log price at time t is x*{t) = pt + l3T{t) + T^/^{t)w{t) where 

Tit) = f v{s)ds = A-i[(l - e-^*)z;(0) + / (1 - e-^'-'~y^)Z{dXy)] 



is referred to as a integrated OU process and models the integrated variance. Quantities of interest 
are often based on the aggregate returns, for s <t, x*{t) — x*{s) which involves 

(3) Tit) - Tis) = A-i[(l - e"^(*-*))i;(s) + / (1 - e-^(*~«))Z(dAj/)] 



where again importantly, w(s) == (e^'^''f(0) + /* e^^(''^^'Z((iAj/)). 

Barndorff Nielsen and Shephard (2001a, Section 5.4.1 and 6.2) show that laws related to the 
random functions 

(4) (^(At),e-^* / c^yZidXy)) 

Jo 

play a key role both in option pricing and likelihood estimation. Specifically option pricing requires 
some type of description of the distribution of 

(5) j\l^e-^^'-y^)ZidXy)^J (1 - e-^)Z(rfy), 

for A = (i — s) > 0. Although the density of (5) is not often known in a nice closed form one can 
apply inversion techniques via its characteristic function or Laplace transform which is described in 
BNS (2001a, 2003). 
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However, as seen in BNS (2001a, 5.4) it is a rather challenging problem to find tractable ap- 
proaches to statistical analysis of likelihood models based on n aggregate returns, Xi — x*{iA) — 
x*{{i — 1) A) over periods of time [{i — 1) A, iA] for i = 1, . . . , n and A > 0.[This framework can be 
extended to intervals of varying lengths say A^]. These models are based on the unobserved actual 
variances Ti — T{iA) — T{{i — 1)A) for z = 1, . . . , n. It is easy to see that one may write 

(6) X, = /xA + Pn + T^/'e^ 

where e^ for i = 1, . . . n are independent standard Normal random variables. Hence it follows that 
conditional on (ti, . . . ,r„) the Xi are independent Normal random variables with unknown mean 
jiA+liTi and variance t^. The major obstacle to tractable analysis of such models is that in general the 
joint distribution of (ri, . . . , r„) is rather complex. BNS (2001a, 5.4) show that statistical inference 
can be done if one were able to sample, efficiently, n iid copies of the pair 

(7) (Z(AA),c-^^ / c>^yZ{d\y)). 

Jo 

The problem is that is not obvious how to deal with the joint distributional behavior of the above 
pair (7). This is in contrast to the option pricing problem which essentially involves the distribution 
of a single random variable. A generic, theoretically all purpose, approach is to use an infinite series 
representation. Several MCMC procedures, based on variations of this idea, have been proposed 
to handle subclasses of these models requiring simulation of points from random processes. See for 
instance, Roberts, Papaspiliopoulos and Dellaportas (2004) and GrifBn and Steel (2005), who use 
compound Poisson process specifications for Z, and the discussion section in Barndorff-Nielsen and 
Shephard (2001a). For approaches to other types of models see for instance Eraker, Johannes, and 
Poison (2003). 

While these methods have their attractive points they do not provide exact solutions for cases 
where Z is an infinite activity process, such as a Gamma process or more generally a Generalized 
Inverse Gamma (GIG) process. Moreover these methods are computationally non-trivial and further 
work needs to be done to assess their accuracy for different processes. Another important point is 
that they cannot be used if one does not have specific knowledge of the Levy density associated with 
Z. This excludes for instance the case where Z is based on a Pareto or LogNormal distribution. 

1.1 Proposal and outline 

This paper focuses on several subtopics related to the issues above. In particular we discuss methods 
that avoid working directly with infinite dimensional components. First, perhaps most remarkably, 
we will show that if one chooses Z to be a Gamma process then one can sample exactly random 
variables based on the pair in (7) and (4). In addition, we will be able to derive the explicit density 
of certain quantities which is also relevant to option pricing. Curiously we will show that the explicit 
densities depend on the dilogarithm function 

r iog(i - u) ^ ^ x^ 

Jo U ^^^ K 

The dilogarithm function is a well-studied special function that arises often in a variety of contexts. 
See for instance Maximon (2003) and Flajolet and Sedgewick (2006). This leads to an explicit 
description of the relevant ti in terms of sums of independent random variables which allows one 
to perform likelihood estimation based on sampling 2n iid random random variables as well as 
the independent random variable f (0). We then easily extend this framework to possibly random 
observation times. An important point is that these results allow one to also use t in other likelihood 
models not discussed in BNS (2001a, 5.4). These facts have not been pointed out in the literature. 
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They are derived from the unique properties of the Gamma/Dirichlet process calculus wherein we 
are able to exploit a, not immediately obvious, connection to Dirichlet Process mean functionals. 
In as much, the seminal work of Cifarelli and Regazzini (1990) and the perfect simulation methods 
discussed in Guglielmi, Holmes and Walker (2002) play a key role. The corresponding OU process 
v{t) is known as OU-F process. This should not be confused with the often discussed F-OU process 
where the BDLP is a compound Poisson process and v{t) has Gamma distributed marginals. Also, 
our results suggest that one could simply use the OU-F as a building block for more intricate models. 

Secondly all the properties that we exploit for the Gamma case do not extend to other OU 
models. However we show that the ability to perfectly sample the marginal distributions of quantities 
relevant to option pricing and likelihood estimation extends to a large class of models where Z is a. 
Generalized Gamma Convolution(GGC). We call these models finite GGC or (FGGC). A highlight 
of this paper related to this class of models is Theorem 6.1. 

Although we shall focus primarily on the BNS OU models, we note that there are many others 
which can be found for instance in Carr and Wu (2004) [see also Carr, Geman, Madan and Yor (2003)] . 
As such we shall employ an analytical technique which leads to an expression of the relevant like- 
lihoods in terms of an n-dimensional Fourier-Cosine integral. This technique is loosely based on 
the ideas in James (2005b). Multidimensional Fourier-Cosine integral appear often in various fields 
including physics. We then focus on ingredients necessary to carry out Monte Carlo procedures 
which are known to be well suited to approximating high-dimensional integrals. More details may 
be obtained from the provided table of contents. 

Remark 1. Throughout, when appropriate, we will be describing the law of a generic positive 
random variable W by its corresponding Levy exponent defined as 

-logS[e-'^^] 

We will use the notation A as an arbitrary positive distance between two points. We shall specify its 
value when necessary, i.e. A ~ t — s, A ~ t and so on. We will also often use the notation a = XA. 

2 Preliminaries 

This paper utilizes results from several linked but not often jointly studied areas. We anticipate that 
the average reader will be familiar with some but not all of the topics. As such we provide some 
details that we shall exploit. The majority of the discussion in sections 2.1-2.2 may be found in 
BNS (2001a, b, 2003). Section 2.3 is again a blend of ideas from several fields. 

2.1 Some more preliminary OU results 

We will describe the distribution of pertinent quantities via their Levy exponents, and discuss the 
basic structure of the likelihood. First note that for any positive g on [0, oo), we may define a random 
variable Z{g) :— J„ g{x)Z{dx). Moreover it is fairly well-known that the Levy exponent of Z{g) 
is given by /„ 'ip{ujg{x))dx. It is clear that all the OU related processes that we encounter are 
representable as some Z{g) where g{x) is readily identified. Using this fact or consulting directly 
BNS(2001a,b, 2003) one has that the Levy exponent of the quantity in (5) is 



(8) / ^{u;{l -u))u-^du 

for A = (t - s) > 0. The Levy exponent of c"^* /J e^yZ{d\y) = q-^^ J^ c^yZ{dXy) for A = t is, 

iP{luu)u~ du 
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If we wish to choose v{t) stationary then the Levy exponent of f (0) must be 

(9) / ip{Luc^'')ds = / il){ijju)u^^du. 

Jo Jo 

2.2 BNS Likelihood model 

The model of Barndorff-Nielsen and Shephard (2001a, section 5.4) translates into a likelihood based 
model as follows. Let Xi for i = 1, . . . ,n denote a sequence of aggregate returns of the log price 
of a stock observed over intervals of length A > 0, described in (6). Suppose additionally the Z 
depends on unknown parameters g. The likelihood of the model depends on unknown parameters 
■d = {iJ,,P,X,g) and as stated before the Xi|?9,T are iid Normal random variables. Ideally one is 
interested in estimating ■(? based on the likelihood 



(10) J-(X|^) = / 



Y[(biX,\^iA + (3n,n) 



fin,. ..,Tn\g,X)dTi,. ..,dTn 



where, setting Ai = {Xi — /iA), and A ~ n ^ 127=i ^i^ 



V Stt 



denotes a Normal density. The quantity /(ti, . . . , r„|£i, A) denotes the joint density of the inte- 
grated volatility based on the intervals [{i — 1)A, iA] for i = 1, . . . , n. Barndorff-Nielsen and Shep- 
hard (2001a) note that the likelihood is intractable and hence makes exact inference difficult. The 
apparent intractability is attributed to the complex nature of /(ri, . . . , t„|£i, A) which is derived from 
a random measure. Specifically, the BNS models complexities arises from the following structure of 
the Ti. From (3) one has for the BNS model 

(11) \n = {l-c-^^)v{{i-l)A)+ f {l-c-^^'^-y^)Z{d\y) 

J(i-1)A 

where importantly for rj ~ e^^ , and 

O, = / e-^'^'^-y'^)Z{dXy), 

J{j-i)^ 



v((i-l)A) ^ e^^^'^^^^[v{0)+Y,]J'-^rjOj]. It is not diffieuh to see that the O^ are iid for j = l,...,n 
but are correlated with corresponding terms 

(1 - c-^^^^-y^)z{d\y) ^ Zj - Oj 

'0-l)A 

where Zj := [Z{XjA) — Z{X{j — 1)A)] = Z{\A). Furthermore O; appears in each Ti for i > I. Hence 
the suggestion by BNS to try to sample the iid pairs in (7). 

Indeed the joint distribution of the (ti, . . . ,r„) is in general complex. However one can easily 
obtain its joint Laplace transform. It is with this fact that we argue that the primary stumbling 
block which currently prevents one from integrating out the infinite-dimensional components in the 
likelihood, is inherent from the Normal distribution of Xiji?, r. Quite simply the Normal assumption 
yields exponential terms of the form 

e-^i/(2n) j-ather than e"'^'^^. 

We will show in the forthcoming sections how to apply a Bessel integral representation, which does 
not depend on the distribution of (Ti,...,r„) to obtain expressions for likelihood based on quite 
general candidates for t. First however we will describe the very remarkable and unique properties 
of the OU-F model in section 3 which does not require this approach. 
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2.3 Some points about GGC and Gamma processes 

Wc will be making extensive use of the basic elements of the theory of Generalized Gamma Con- 
volutions (GGC) which can be found in Bondesson (1979, 1992) and Thorin (1977). GGC are a 
sub-class of infinitely divisible random variables. A nice point is that they all have the important 
self-decomposability property. This has an interesting consequence since it is well known that v{t) is a 

stationary OU process if and only if v{0) = v{t) is self-decomposable. See for instance Wolfe (1992), 
Jurek and Vervaat (1983), Sato (1999), Jeanblanc, Pitman and Yor (2002) or BNS (2001a, Theorem 
1) for a more precise statement. That is, there is a large subclass of OU models which all have GGC 
laws. Some important examples of GGC random variables and corresponding processes are GIG 
laws. Stable laws of index < a < 1, and of course Gamma random variables. 

Important, from our point of view, is that a random variable is a GGC if and only if its Levy 
exponent is expressible as 

/■OO 

(12) / de{ujx)iy{dx) 

Jo 

for some arbitrary sigma-finite measure satisfying appropriate conditions so that (12) is finite and 
where 

/•OO 

(13) dg{u;) = 9\og{l+uj)= {1 - c-'^'')9s-'^e~-''ds. 

Jo 

corresponding to the Levy exponent of a Gamma random variable with shape parameter 9. That is 
to say dg is a special case of ip and moreover the Levy density of a corresponding Gamma process 
is given by 

peids) = 9s^^c^''ds for s > 

It then follows that the Levy density of a GGC is given by 

/•OO 

(14) Os-^ / c-'''v{dr). 

Jo 

As a consequence, if we denote a Gamma process on a Polish space ^ with sigma finite shape 
measure 9v* as Ggi,*, then (12) is significant as it coincides with the Levy exponent of an arbitrary 
Gamma process mean functional say Ggiy* (g) — J^ g{x)Ge^* {dx) where by a change of variable, 
R — g{X), one can write equivalently in distribution as /^ rG0^{dr). 

Now we point to a key fact that has not been exploited much in the literature. First throughout 
this paper let Tg denote a Gamma random variable with shape 9 and scale 1. Denote the density of 
a Gamma random variable with shape 9 and scale 5 > as 

When & = 1, we simply write '^0(1/). Let (Ji) denote the jump points of a Gamma process and let 
(Zi) denote the points of a Poisson random measure whose laws are determined by v, which are 
independent of (Ji). It is well known that one can write Gouidx) — X]"=i Ji^Ziidx). Furthermore, it 
follows that if y is a GGC random variable then one can always write 

Y^Ge.49)^TeMe, 

where Mg^, — '^i^i{Ji/Tg)Zi, is a random variable independent oiTg. The independence property is 
due to the known fact that the sequence (Ji/Tg) of probabilities is independent of Tg which may be 
written as Tg — X]i=i Ji- This property uniquely characterizes a Gamma process and has nothing 
do with whether or not v is finite or more generally sigma finite. The sequence {Ji/Tg) is known to 
have the Poisson-Dirichlet law. 
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Hence when ly :— H is a finite measure, which we will take without loss of generality to be a 
probability measure, a Dirichlet Process with shape 9H, having total mass 9H[^) = 6, is defined 
by the representations 

where importantly To = GeHijX') is independent of Pen- Setting ^ ~ [0, oo) one has that 

f-OO 

Men = / xPenidx) 
Jo 

is a Dirichlet Process mean functional which again is independent of Tg. This independent property 
naturally comes from the finite dimensional Beta-Gamma calculus based on the classic result of 
Lukacs (1955), which we shall also use. That is, if Tg. for i = 1, . . . ,n are independent Gamma 
random variables with shape 9i then the sum, X]i=i ^^i ~ Tg^,, where 9* = X^^i ^i ^^^ moreover is 
independent of the vector of probabilities (Tg./Tg*) which has the Dirichlet distribution with density 

n 

D{pi,...,Pn) OC ]Jpf"^ 
1=1 

where the Y!^^=iPi — 1- ^'^ ^'^^^ denote the fact that a random vector has Dirichlet law of this type 
by writing DiRiCHLET„(0i, . . . , 6'„). Similarly denote a two parameter beta law as Beta(0i, 6*2). 

2.3.1 Connection to Cifarelli and Regazzini distribution theory 

Because of these observations we are able to exploit the works of Cifarelli and Regazzini (1990) and 
those of subsequent authors to obtain expressions for the marginal densities of relevant components 
of large class of models which we call OU-FGGC. The FGGC are models with Levy density defined 
by (14) with v := H. This is relevant to both option pricing and likelihood estimation. We should 
add that many of these properties will extend to more general moving average models where Z is 
an FGGC BDLP. 

The study of properties of Dirichlet process mean functional has been a major area of interest 
in Bayesian Nonparametrics. This line of work was initiated by the paper of Cifarelli and Regazz- 
ini (1990). One of their important contributions was to obtain explicit expressions for densities of 
mean functionals MgH- Let /mo^ denote the density of Mgn- Set H{x) = L H{du). Then from 
Cifarelh and Regazzini (1990) or Cifarelh and Mehlh (2000) one has for = 1 

(15) Jmh [x) = - sin(^i7(x))e- /o°° iog(l*--l)H(dt) 

TT 

and when 9 > 1, 

(16) fMeH{^) = ^—-^ I (a;-w)''"^-sin(7r6'iJ(M))e-''''o°°'°s(l*-"l)^(''*)du 

T^ Ja 71" 

One can also obtain an expression for the cdf of Mgn that holds for all > 0, we do not list that 
here. 

2.3.2 Perfect simulation of Mgn 

It is evident that given the form of the density in (15) one can in principle use some sort of rejection 
sampling procedure to obtain realizations of Mh- With a bit more care one can devise an efficient 
method to sample Mgn for 9 > 1 using the density in (16). Importantly, as pointed out by Hjort 
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and Ongaro (2005), when 6 ^ m, where tti = 2, 3, . . ., is an integer one can use (15) to sample MmH 
based on the following fact, 

m 

1=1 

where (Mi.i) are iid with common distribution equivalent to Mh given by (15). Moreover (Pi , . . . , Pn) 
is independent of (Mi^j) and is n-dimensional Dirichlet„(1, . . . , 1). This can be seen as a simple con- 
sequence of the infinite divisibility oiTmMmH, where , as a consequence, TmMm.H — J2^i Ti.iMi^i, 
and applying the Beta Gamma calculus. That is further writing Pi ~ Ti^i/Tm, where Tm — X^i^i ^i j 
is independent of (Pi). What is important is that these methods do not rely on the more compu- 
tationally burdensome, and otherwise approximate, series methods. There is however yet another 
approach which will allow one to easily perfectly sample Mgu for all > 0. 

Recently, in the case where Mqh is almost sure bounded, Guglielmi, Holmes and Walker (2002) 
devise a very simple and efficient method to obtain perfect samples from the distribution of Mqh that 
works for all 6 > 0. We recount the basic elements of that algorithm. First note that < a < Mqh < b 
if and only if the support of H is [a, 6]. As explained in Guglielmi, Holmes and Walker (2002), 
following the procedure of Propp and Wilson (1996), one can design an upper and lower dominating 
chain starting at some time — A^ in the past up to time 0. The upper chain, say uMgn, is started at 
uMgH,-N = b, and the lower chain, IMgn, is started at IMgn.-N ~ ei. One runs the Markov chains 
for each n based on the equations, 

(17) uMeH,n+l = Bn.eXn + (1 - Bn.e)uMeH,n 

and 



(18) IMgn^n+l — Bn.eXn + (1 — Bn.e)lMgH,n 

where the chains are coupled using the same random independent pairs (P^.e, Xn) where for each n, 
Bn^g has a Beta(l, 9) distribution and Xn has distribution H. The chains are said to coalesce when 
D ~ \uMgH,n ~ IMgH.nl < £ for some small e. Notice importantly that this method only requires 
knowledge of the distribution H. 

Remark 2. Vershik, Tsilevich and Yor (2004) and James (2005a) are two examples of applica- 
tions that directly exploit the independence property exhibited at the level of the Gamma/Dirichlet 
process. See also Diaconis and Kemperman (1996) and Diaconis and Freedman (1999) for more 
interesting facts. 

Remark 3. More discussion on the merits of self-decomposability as it relates to financial 
applications can be found in Carr, Geman Madan and Yor (2005). 

3 Laws and Likelihoods for the OU-F model 

For 9 > 0, define a OU-F process by setting Z = Gg, where Gg denote a homogeneous Gamma 
process on [0, oo), i.e. v{dx) = dx for x e [0, oo) with law specified by its Levy exponent dg{uj) given 
in (13). Letting vg{t) denote the stationary OU-F it follows that its Levy exponent is 

(19) / dg{uju)u-^du = 9 (l-e-'^y)y-^Ei(y)dy = -0Li2(-uj) 
Jo Jo 

where Pi(y) = / c^^u^^du — J. e^^^u^^du is Euler's exponential integral. That is to say the 
Levy density of vg{0) is Pveidy) = 6y^^Ei{y)dy. 
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Remark 4. In addition to obtaining the form of the Levy density, BNS (2003, p. 283) note that 
the Levy exponent of a OU-F can be expressed as. 



^(-1)^^ forO<tj< 1 



but they don't equate this with the dilogarithm function. 

The previous discussion indicates that one can implement both option pricing and likelihood 
analysis if one can sample the special case of (7) given by 

.A 

(Ge(AA),e-^^ / c^«Ge(dAy)). 
Jo 

The Levy exponent of the second term is given by 





/ dg{uju)u ^du — 1 dga{'-^u)Fa{du) ~ —9[Li2{—u}) - 
Je-" Je-" 


- Lz2(-a;e-»))] 


where 






(20) 


jQ-a au a 





is a cdf for e " < y < 1. However due to the fact G0(AA) — GgaPa ([e ",!])= Tga this is equivalent 
to sampling the pair 

{Tga, / xPgaF^idx)) 

Jc-" 

where for e~° < y < 1 

p ,j X GgapAdy) 
PeaFAdy) == ^^ 

J- 9a 

is a Dirichlet process random probability measure with shape parameter OaFa- 

Remark 5. We shall use the notation Mga rather than the perhaps more accurate MgaF^ where 
it is understood that Fa is defined in (20) 

We discuss some of the implications of these facts in the next two propositions. 

Proposition 3.1 For each fixed A > 0, and a = AA set Yga '■— e " /g e^^Gg{dXy), where Gg is a 

homogeneous Gamma process. It follows that Gg{a) ~ J^ Gg(d\y) — Tga- Additionally, the following 
distributional properties hold. 

(i) Let Mga '.— L xPgaFa{dx) denote a Dirichlet process mean functional based on the shape 
parameter OaFa. Then for each fixed A, one has the coordinate-wise equivalence in joint dis- 
tribution, 

{Gg{a),Yga)^{Tga,TgaMga) 

where Mga is independent of Tga. Furthermore e^"" < Mga < 1 almost surely. 

(ll) /„*(! - e-^(^-^))G9(dAy) == Gg{a) - Yga = Tga[l - Mga] 
(ill) {Gg{a) - Yga, Yga) = {Tga[l - Mga],TgaMga) 
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Proof. The result is already established by our construction and appealing to the unique 
independence property of the Gamma/Dirichlet process. However since the joint equivalence in 
statement (i) is the key factor separating the OU-F from other OU processes, hence quite delicate, we 
will check it via joint Laplace transforms. Evaluating the joint Laplace transform of the {G0{a), Iga) 
at points {llIi,ll!2), it is easily seen that the joint Levy exponent is 

1 

doai^l +UJ2u)Fa{du). 

Now being careful to use only the independence property of Tga and Mga and the fact that Moa 
is a Dirichlet process mean functional we proceed as follows. Write LOiTga + ^2TBa^ea = Tga[^^\ + 
W2Mea] := W. Furthermore note the loi + uj2Mga = /q (wi + uj2x)PgaF^{dx) = PgaFaid), for g{x) — 
LOi + L02X. Now by independence of Tga and Mga the joint Laplace transform, taking expectation 
with respect to the Gamma law first is, 

Eie-'^] - E[(l + ^1 + u:2uMga)-'^] = E[(l + PeaFaig))'""] 

Now appealing to the well-known identity of Cifarelli and Regazzini (1990) it follows that 



E[(l+PeaFa(g)) 



— Wai 



= in, C 



-GeaFcia)] 



which is the desired result. The above argument indeed establishes the proof but the very special 
nature of the result perhaps will not be fully clear until one reads section 3.8. □ 

The next result describes the distribution of vg{Q) in the stationary case. 

Proposition 3.2 Let vg{Q) have distribution described by the Levy exponent (19). Let Gov denote 
a (non-finite) Gamma process on [0, 1] with v{du) — u~^du where L v{du) — 00. Then vg{Q) is a 
generalized Gamma convolution (GGC) such that 

/■I 
ve(Q)^ / xGeu{dx)^TeMg, 
Jo 

where Mg — Mg^, is independent of Tg but is not a Dirichlet process mean functional. Furthermore, 
for each fixed 0, the distribution of Mg is characterized by its generalized Cauchy-Stieltjes transform, 



Remark 6. It is quite possible to obtain an explicit form of the density of vg{0) by using 
standard inversion results for characteristic functions and noting the relationship of the complex 
valued dilogarithm function to the Inverse Tangent Integral, 

„ , , f^ arctan(u) , 

Ti2{y) = / —du, 

Jo " 

which is the imaginary part of the complex valued dilogarithm function, and Clausen's Function. 
For more details see Maximon (2003). 
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Recapping, Proposition 3.1 shows that the distribution of (Ge(a),l0a) is determined by the 
distribution of the independent random variables (Tga, M^a)- Among OU processes discussed here, 
the independence property is unique to OU-F processes. Additionally, as we shall see this pair may be 
sampled exactly due to the fact that MgaPa is & Dirichlet process mean functional. On the other hand 
Proposition 3.2 shows that although vg{0) = TgMg is a GGC, the results for the Dirichlet process 
do not apply to Mg and we otherwise do not have a tractable expression for the explicit density of 
vg{Q). However, we do believe that a careful use of the relationships mentioned in Remark 6 will lead 
to an explicit form. The next proposition, using the work of CifarcUi and Rcgazzini (1990), provides 
more details for the distribution of Mga and shows also that one can use the Dirichlet process results 
to obtain a good approximate for the distribution of we(0). 

Proposition 3.3 For each < a ~ AA < oo and 9 > 0, let Yga = e^^^ L e^yGg{d\y) denote an 
infinitely divisible random variable with Levy exponent, 

1 .1 

dg{ijju)u du — I dga{uJu)Fa{du) ~ —6[Li2{—oj) — Li2{~oje "))] 

where Fa{du) = a^^u~^du is the density of a random variable taking its values in the interval 
[e~", 1]. Then the following results hold 

(i) Yga = TgaMga, where Mga — /g-a xPgaFa{d,x) IS a Dirichlet process mean functional. 

(ii) The Levy density of Yga is pga{dy) =^ 6y^^[Ei{y) — Ei{ye"')]dy. Hence the cumulants of Yga 
are for each integer j , 

9 r y^-'[E,{y) - E^{ye'')\dy - 0^(1 - g-"^) 
Jo J 

(Hi) When 9a = 1, the density of Mi is given by 



(21) - sin 

TT 



-7rlog(a:;) 



^^ll^logix)]^l^f^^^[Li2ix) + Li2(^)] 



for e "■ < X < 1. 
(iv) When 9a ~ \, the density ofVa := — log(Mi)/a is given by 

(22) isin(7rz;)e-[-+-'le5^e^[^^^(e""")+^^^(e""'""')l, 

TT 

/or < V < 1. 
(v) When 9a > 1, the density of Mga/ a is given by 

(23) ^^^ /' (x-e-"")'""'sin(7r0a«)e-^"["+"'le^e-^[^'^(e"°'')+^^^(e" 



- log x/a 

forO<v <1. 

(vi) Lf 9a — m, where m — 2,?!^ ... is an integer, then Mm ~ X]i=i WiMi^i, where (Mi.i) arc iid 

with density (21) and independent of (Mi^i), {Wi — Ti^i/^-^ 2^i,i); where , Ti^i = Ti are iid, 
is a DlRICHLET„(l, . . . , 1) n-dimensional vector. 

(vii) Yga converges in distribution to vg{0) as e^" -^ 0. 
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Proof. Most of the results are immediate from our previous discussioms. The forms of the den- 
sity arises from apphcation of CifareUi and Regazzini (1990) which amounts to explicitly calculating 
/e-a log(|t — x\)Fa{dt) expressed in terms of the dilogarithm function. D 

The last result in this section gives a completely tractable description of the conditional distri- 
bution of the log asset price at time t given information up to time s. This type of result is pertinent 
to option pricing as discussed in BNS(2001a, 6.2) and Nicolato and Vcrnados (2003). 

Proposition 3.4 Let Xg{t) be defined as in (1) with Z = Gg- Additionally for Q < a < t, set 
A = (t — s) and define ft,(A, s) = (1 — e^^ )v0{s) and jj,'* = /iA + Xg(s) + (3h{A,s). Then the 
conditional density of Xg(t)\xg{s),vg{s) is given by 



(/-(a;]^* + (3y, /i(A, s) + y)qea{y)dy 
where qea{y) = /e-a %a(y|(l - u))/m(,„ (w)dw. When 0a = 1, 

gi(y)- / ^i{y\{i^e--))fvMdv 
Jo 

where fv^ is the density of Va given in (22). O 

3.1 Perfect simulation of Moa 

Due to the fact that the dilogarithm function, Li2{x), is a well- understood special function, which 
is available in many computational packages, it is evident that the densities in (21) and (22) can be 
exactly sampled using a rejection procedure. Again based on the discussion in section 2.3.2 Statement 
(vi) of Proposition 3.3 shows that one can use this fact to easily obtain samples of Mm, and hence 
Ym, for any integer to. With a bit more care one can devise an efficient method to sample Mga for 
6a > 1 using the density in (23). One can also use the perfect sampling method described in 2.3.2 for 
all 6a, based on uMga-N = 1 and lMga,-N — c^", Bnga is Beta {1,6a) and X„ has distribution 

Fa 

3.2 BNS OU-r likelihood inference 

The results in the previous section now give the ingredients to perform likelihood based statistical 
inference via simple exact sampling. Here we describe a bit more about the distribution of r^ in the 
OU-P case and then extend the discussion to randomly sampled times. 

Proposition 3.5 Define for A > 0, a = AA and i = 1, . . . ,n Tgj :— Tg{iA) — Tg{{i — 1)A), by 
setting Z = Gg in (3). Furthermore, let ri = e'^*^ for i = 1, . . . ,n, with tq = 1. Then it follows 
that, for i ~ 1, . . . , n, 

(24) Are., = (1 - e-^'^)vg{{i - 1)A) +T,[1- M,] 

with, 

i-l 

vg{{^-l)A) =. e-^(-i)^[z;,(0) +5]r,r,M,.] 
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where {Ti,Mi) are iid pairs independent of vg(0). Additionally, for each fixed i, Ti and Mi are 

independent with distributions specified by Ti — Tga and Mi = Mga- This implies that likelihood 
inference for the model (10) may be obtained from the joint distribution of {Xi,Ti, Mi,vg(0)) given 
by 



(25) 



Y[ (j){X,\fiA + PTe^,,T0^,)%a{U)fM{vi) 



.1=1 



fve{Q){w) 



where Tg^i is expressed as in (24), with Ti = ti,Mi ~ Vi, and vg{0) ~ w. O 

A Bayesian procedure, which involves placing a prior on ?9 == (/i,/3, A, g), is quite natural and oth- 
erwise proceeds by standard arguments, in this setting. That is letting 7r(??) denote a prior joint 
density it follows that a posterior distribution of x?|X is determined by a posterior distribution of 
■d, {Ti, Mi), we(0)|X which is proportional to 



7T{d) 



Y\_<l>{Xi\lJ,A + PTe^i,Tgs)%a{ti)fM{Vz 



fve{0)iw) 



Remark 7. The likelihood in (10) for the OU-F case obviously is obtained by integrating out 
the pertinent independent quantities in (25). Due to the Gamma distributions, the answer could be 
expressed in terms of integrals with respect to modified Bessel functions. Or otherwise a subclass of 
Generalized Inverse Gaussian(GIG) random variables. 

Remark 8. Note that in practice we can approximate a draw from the distribution of vg{0) 
by using instead Ygs for e^ small. Otherwise, if strict stationarity vg(t)is not a concern, one can 
certainly use any positive distribution for vg{0). 



3.3 The likelihood via a connection to Variance Gamma processes 

Recall in the stationary case that according to Proposition 3.2. vg{0) — TgMg, where Mg is not 
a Dirichlet process mean functional. However this point allows one to write rg and (re,!, . . . ,Tg^n) 
in terms of a product of a Gamma random variable and another independent random variable. 
Specifically, for a = A A, one may write 



Te,i — Tg(^i+na)Si 



where for i = I, . . . ,n 

A5, = (l-e-°)e-"('-^) 



Tg 



T 



e{l+na) 



M8+J2 



T, 



— Tg(l+na) 



r,M, 



T 



T 



Mi 



e{l+na) 



The vector S = (^i, 
may also write 



, Sn) is independent of ^^(i+na) which can be written as Tg + X]r=i -^*- ^^ 



AS*, = (1 - e-")e- 



-a(i-l) 



Pn+lMg + ^PjrjM-j 
J = l 



P^[l-M,]. 



where Pn+i = 1 — X]?=i Pj^ ^^^ (-^i' • • • ' Pn+i) is DiRiCHLET„+i(6'a, . . . , 6a, 6) independent of all 
other random variables. Recall now that a GIGiy, 5, 7) random variable has density given by 



g{x\v,5,"f) 



2K,{djy 



.-l^-i(5^.-H7^x) fo, ^ > 
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where K^, is a modified Besscl function. Recall also that K^{x) = K-^{x). 

Additionally we will exploit the following nice feature of K^{x). Suppose that for a to = 0, 1, 2 . . ., 
the |:^| = TO, + 1/2, where \v\ denotes absolute value, then we can use the fact that 



(26) 



K„ 



.+ 1/2 



(x) 



E 

fc=0 



[m + fc)! 
k\{m - k)\2*'^ 



See for instance Pitman (1999, eq. (40)) for a probabilistic interpretation of (26). 
This facts leads to the following description of the likelihood. 

Theorem 3.1 The observations according to (6) can be represented as X^ = /iA + PTgiij^_na\Si + 



nl/2 



Te(l+na)S^Y'e,, m the OU-T case. Setting 7^ ^ [2 + P^Y.i^^S^] and 5^ ^ E"=i ^'/(25'j), k 
?(1 + na) and v = k — n/2, and a ~ AA. The following results hold. 

(i) The likelihood in (10) can be written as, 



.if(X|i?) 



e"'^'^Ei, 



2K,{6-i) 



n 



(ii) If 9 and a are chosen such \v\ ~ m + 1/2, for to, = 0, 1, 2, . . ., then 



^(X|^) = e^^f' f2 ...^""^fj.L lE^ 



k=0 



fc!(TO-fc)!2'= 



(7/<5)"r(«) n V2^^« 



As a special case \v\ = m + 1/2 for all n, if 9a = 1/2 and 9 = m + 1/2. 
(Hi) If additionally m — 0, that is 9 — 1/2 and a = 1, then 



^(X|z9) == e"^%^ 



-5.27- 



n 



r(K) fj^ V2^ 



5, 



-1/2 



In all cases the distribution of (6*1, . . . , S'„) is completely determined by the 2n + 2 independent 
random variables with joint density [[\i^'^9a{ti)fMea{'>^i)]'^9{t)fM (^)^ 

The next result in effect serves to make clear Theorem 3.1 but also highlights the possibility, from 
a practical point of view, for more data augmentation procedures 

Proposition 3.6 Consider the setup and notation in Theorem 3.1. Additionally define 01 — 0^ X]?=i ^i 
Then it is clear that 






y-te-^(^^«"^+««)^.(y)dy, 



which leads to another expression of the likelihood ^(X|i9). Thus statistical inference may be based 
on simulation from the joint density 



^M 



^%a[U)fMeAv^) 



Mt)f^^{w). 



Based on this fact one has that if a random variable V has the density '^k relative to the representa- 
tion o/^(X|'i?) then the posterior distribution of V\S, X, d is GIGiy, 5, 7) with parameters specified 
by Theorem 3.1. and Proposition 3.5 
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Remark 9. One notes that the expressions in statements (ii) and (iii) of Theorem 3.1 arc quite 
manageable. Here one is perhaps taking the view that 6 and AA are chosen to case computations. 
However note that in statement (ii) that to, whose parameter space is {0, 1, . . . , } becomes a viable 
and flexible parameter of interest from a modelling point of view. The expression in statement (i) is 
also quite amenable to Monte-Carlo estimation approaches. 

Remark 10. By Variance Gamma processes we are loosely referring to the work of Madan, Carr 
and Chang (1998), see also Carr, Geman, Madan, and Yor (2003). It is evident that all OU-GGC 
models exhibit similar properties. That is if the BDLP Z is a GGC then analogues of Theorem 3.1 
and Proposition 3.6 have exactly the same form. However, in contrast to the OU-F case, one still 
does not have an obvious way to sample from the distribution of S. 

3.4 Bayesian estimation and related comments 

We have shown that the distribution of (rg i, . . . , rg „) is determined by 2n + 2 independent random 
variables whose distributions can be perfectly sampled or in the case of we(G) approximated with 
arbitrary accuracy. We also note that the explicit densities that we have given for Mga definitely 
have practical utility, whereby rejection methods can be used. We also believe they are interesting 
from a mathematical point of view as they may have connections to application in physics or analytic 
combinatorics. These are places where the dilogarithm function appears often. However, in terms 
of practical simplicity it is perhaps easier to use the perfect simulation schemes which work for 
all values of 9a and only require simulation from beta random variables and the distribution Fa- 
Also, in regards to vg{0), we note again that in the case of not strictly stationary OU-F models, 
we may choose vg{0) to have any distribution. However Theorem 3.1 suggests there are some quite 
interesting simplifications that occur if we choose ^^(G) = TqW , where W denotes a random variable 
independent of Tg. We note again that all GGC random variables have this form including the class 
of GIG models. 

Armed with the information that we have provided one can construct a variety of efficient sim- 
ulation based techniques. Here we briefly highlight the Bayesian approach. Primarily this is due to 
the fact that a Bayesian approach is essentially an approach involving integration and hence is a 
quite natural for Monte-Carlo based estimation. It is in many respects quite similar to Bootstrap 
techniques. We now mention some well known points about Bayesian estimation. Suppose that ■n{d) 
is a prior distribution of the unknown parameters. Then, as is well known, the fundamental object 
of interest is to obtain the posterior distribution of 0|X, which is given by 

7r(t?|X) (X 7r(i9)^(X|i9) 

Estimation of some parameter h{{)) can then be cast in terms of integration, 

^/,qV«Ag 2K^{Si) -r-rn 1 0-1/2 

"■we (7/5)T(k) 11*=i ^n2^^i 



(27) E[h{d)\X\ = / h{u)T:{du\X) 

Je 



where the denominator should be understood as. 



riAd 2K,4S-t) Yfn 1 c-1/2 



^(X) = / 
Je 



e 



nAl3 2ir^((57) -A- _J_ 1/2 

(7/<5)^F(«) IJ V2^ ^ 



7r(di?). 



For instance, the posterior probability that ?? is in some region B can be evaluated by choosing 
h(x) = I{x G B}. Since Bcssel functions, such as Ku{x), are available in standard mathematical 
computer packages, one can just draw from the joint distribution of (??, S), which is readily available 
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from our results. That is for I = 1, . . . , _B draw iid random vectors (?9;, 5*1.;, . . . , Sn.i) then (27) is 
approximated by 



(28) 



1.1=1^ 1^,/5,YT(kA i.i.i=l ,/2^'^iJ 



The nice feature of basic iid Monte-Carlo type estimator like (28) is that accuracy issues are well- 
understood and are less dependent on the sample size. Here accuracy increases as B increases. 

One can of course develop more sophisticated importance sampling and MCMC methods based 
on well-known ideas. These may involve sampling from the posterior distributions. For instance, our 
results show that the posterior distribution of i^jX can be obtained by working with the posterior 
distributions of t^lX, S, V and S, V^|X, ■& where for instance V\S, X, d has a GIG{v, (5, 7) distribution. 
All other conditionals can be easily deduced by various augmentations of the expressions given in 
Theorem 3.1 and Proposition 3.6. 

3.5 OU-F processes with possibly random scale parameter 

Up to this point we have assumed that Gg was a homogeneous Gamma process with scale parameter 
equal to 1. This was done mainly for notational convenience. However, it follows from our analysis 
that the introduction of a scale parameter say C, can be used as a powerful modeling tool. Naturally a 
scale parameter can just be introduced by replacing Gg with QGg throughout. However an important 
fact is that if we use QGg, the vector S described in section 3.3 still does not depend on C,. This 
means that one can now write 



Xi — /iA + KTg{l+na)Si + [C,Tg{^i^na)Si 



1/2 



e,; 



Note that if Q is fixed then all our results carry over without change. This means extending the 
model to the case where Q is random is straightforward. The main feature being that we would now 
be working with a Gamma scale mixture, based on C2^ei(i+na)j which can be used to introduce more 
distributional modeling flexibility. 

3.6 Likelihoods for Superpositioned OU-F 

BNS (2001a, p. 178) propose the idea of superpositions of independent OU processes to alter the 
auto-correlation structure. Here, letting p denote a positive integer, and {wi, . . . ,Wp) a possibly 
unknown vector of positive terms summing to l,we discuss briefly a generalization of Theorem 3.1 to 
the case where one starts with a superposition process v{t\p) = X]?=i ^j^Sj (*) where for j = 1, . . . ,p, 
vg-{t) are independent OU-F processes which are based on parameters {\j,0j), in place of {X,9). 
Obviously the distributional results we have developed apply to each of the independent components. 
One uses for instance Oj = AjA and 0jaj in place of a and Oa. 

Let T{t\p) = 'l2'j=iWjTg-{t), denote the integrated volatility where each Tg^{t) ~ J^vg-{s)ds. 
Additionally the analog of (rgi, . . . , Tg^n) is r^ :— T(iA) — T((i — 1)A). Then by similar arguments 
to the previous section one can write for ^„ — X]?=i ^j(l + nXjA), 

where Si^p := Ti/T^^ has an obvious description by applying our previous results to each component 
Tg- , and the vector (S'l^p, . . . Sn.p) is independent of T^^. 

Proposition 3.7 Let Xi = /iA + PT^^Si^p + [T^^Si^p] e^. with terms defined in this section. Set 
72 = [2 + /32 J2"=i S,,p] and S^ = E]Ci AJ/i2S,,p), k - Cn = Ej=i ^j(l + ^XjA) and v = k - n/2. 
Let -dp denote the enlarged parameter space containing unknown quantities such as {wi, . . . ,Wp), 
then the likelihood ^(X|?9p) has the same form as the likelihood in Theorem 3.1 with appropriate 
substitutions of the above parameters and (S'l^p, . . . , Sn,p) in place of {Si, . . . , S'„). In particular, 
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ft) if ELi ^j Aj] A ^ 1/2 and Yf,=x Oj=m+ 1/2 for m == 0, 1, 2, 



then 



■^(xI^p! 



_ nA/3 



E 

fe=0 



(to + fc)! 
A:!(m-fc)!2fc 



E 



— fe _i n 



- (7/<5)"r(..) 



1 

n4 



2tt '^^ 



1/2 



T/iJs expression holds more generally for v = m + 1/2 or v = —to, — 1/2. 
(a) If additionally m = 0, that is Y^^=i ^j = 1/^ '^"■'^ [X]?=i ^j'^j]^ = V^; the 



^(X|z?p) = e"'^''E 



-57 



27^1 

r(«) 



n 



27r^^'^ 



-1/2 



D 



Remark 1 1 . Note that superpositioning allows more flexibility in terms of the parameter values 
for the constraints J2^=i^j = rn + 1/2 and E^=i^i-^j]^ = 1/2- But otherwise preserves the 
simplicity of the likelihood as seen in (i) and (ii) of Proposition 3.7. 



3.7 Randomly sampled times 

From a practical point of view it may be desirable to sample at uneven or random intervals. See 
for instance Ait-Sahalia and Mykland (2003, 2004). The next result shows that the independence 
structure still holds (conditionally) but that the individual terms are not identically distributed. 

Proposition 3.8 Let = 70 < 71 < 72 < ... < 7n denote n random times and define Ai := 
ji — 7i-i. Define Tgj := T0{ji) — Tg{'yi-i), and ri — e ''' for i — 1, . . . ,n, with ro = 0. Then it follows 
that, conditional on (Ai, . . . , A„), for i = 1, . . . , n, 

i-l 

\Te,^ = (1 - e-^^^)e-^^^-' [vg{0) + ^ r,T,-M,] +T,[1- M,] 
where {Ti,Mi) are conditionally independent pairs independent ofvg{0). Additionally, for each fixed 



i, Ti and Mi are independent with distributions specified by Ti = Tgx^. and Mi 



M, 



{e\Ai)F^^ 



If 



the Ai for i — 1, . . . , n are independent then the unconditional distribution of the pairs (Ti, Mi) are 
independent. 



3.7.1 Time changed Integrated OU-F processes 

Notice that the previous proposition places minimal constraints on the possibly random times (7^). 
Naturally if one can easily sample (Ai, . . . , A„), then this would lead to models which are amenable 
to likelihood estimation. These observations lead us to introduce briefly a class of time changed 
integrated OU processes defined as 

rZ(t) nZit) 

(29) re{Z{t)) = / v{s)ds = X-^[{1 - e-^^(*))«9(0) + / (1 - c-^(^(*)-^))G9(dAy)] 

Jo Jo 

where Z is any subordinator independent of Gg. The next result shows how this model is represented 
by Proposition 3.8. 

Proposition 3.9 Consider Tg{Z{t)) defined as in (29). Fori ~ 1, . . . ,n, defincTgiz '■— Tg{Z{iA)) — 
Tg{Z{i — 1)A)). Then it follows thatrg^i^z is equivalent to a specific Tgi in Proposition 3.8 by setting 

7j = Z{iA). Furthermore A^ == Z{iA) - Z{{i - 1)A) = Z(A) are iid.U 
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Remark 12. The time changed process (29) represents an extremely rich class of models which 
adds a great deal of distributional flexibility to the OU-F models. As seen from Proposition 3.9 
likelihood analysis for such models is again easily accomplished. In that case there may be additional 
unknown parameters associated with Z. For instance selecting Z such that 

Corresponds to the case where Z(A) is an Inverse Gaussian random variable. 

Remark 13. One may also replace Z{t) in (29) with any tractable increasing process. For 
instance one may choose t* (i) to be an integrated OU-F process independent of Tg 



Remark 14. Leverage type models discussed in BNS pose no extra difficulties. In the simplest 



1 /2 

likelihood setting, this translates into replacing Xi = /iA + /St^ + t^ e^ described in (6), with 



Where T^ — Tga and r^ is otherwise related to T^ by the representation given in Proposition 3.5. v 
is a real- valued unknown quantity. 

Remark 15. We can extend the OU-F processes based on the homogeneous process Gg to one 
based on an inhomogeneous Gamma process Gg^, where v is an appropriately defined sigma-finite 
measure. That is the Levy exponent for any positive function g of Gg^{g) = L g{x)Gg^{dx) is given 
by /p dg(ujg(x))v{dx) . The volatility process is then defined by 

vg,.{t) = e-^*«(0) +e-^* / c^yGgxu{dy) 
Jo 

The process vg^{t) is stationary only in the homogeneous case. However the independence properties 
that we exploited still hold and one has fairly obvious generalizations of the results we have presented. 
An advantage is that this is another way to increase distributional flexibility. 

3.8 The special nature of the OU-F process as an SV model 

It is important to note that this independence phenomena, exhibited in Proposition 3.1, which allows 
one to easily describe the joint structure of (ti, . . . , r„) for a potential SV model is not only due to 
the usage of a Gamma process Gg. That is to say it will not necessarily be true for non-OU models 
based on Gg . To see this deflnc a moving average process of the type 

t 
(t-x)e- (*-=") G(,(dx) 

It is not difficult to see that the analog of (7) amounts to {J^ c^yGg{dy), J^ ye^yGg{dy). To see the 
problem first set Ha to be uniform [0, a], and gi{y) = e^^, and .92(2/) — y^^- Then it clear that the 
pair above are equivalent in distribution to the pair 

(30) {TgaPgaH^{gi),TgaPgaHa{92)) 

The good point about this representation is that the marginal distributional results for Dirichlet 
process mean functionals apply. This means, for instance, that basically all Levy moving average 
processes that are driven by a Z which is an FGGC have the property that any calculation involving 
a one-dimensional random variable can be calculated using the marginal distributional results for 
Dirichlet process mean functionals. This has an immediate consequence for option pricing formula 
based on such models. 



20 OU-Gamma 

However it is quite clear from (30) that one can negotiate the dependence structure in a manner 
similar to Proposition 3.1, if and only if PgaHaiSi) can be expressed as a function of PeaHaiQ^)^ 
which is not true for this example. This is also why in (7) the OU-FGGC models we shall discuss 
do not have the structure exhibited in Proposition 3.1. In other words Z{a) in that expression has 
to have a Gamma distribution. Or more generally expressible as Tga and a function of the other 
coordinate. Of course the OU-F is not the only Gamma driven SV model that has the ability to be 
exactly sampled as we did in this section. Another example is the Dykstra and Laud (1981) type 
model, see also James (2005b, p. 1784, eq. (29)), which takes the simple form 

rt 

{t-x)Ge{dx). 



In this case the analogue of (7) amounts to {Gg{a),L yGg{dy)). 

4 General Likelihoods 

We now proceed to show how one may perform likelihood analysis for more general (ti, . . . , r„) 

4.1 Fourier-Cosine integral representation of the likelihood 

In order to calculate (10) we use the classical Fourier-Cosine integral 

1 /"OO 2 1 A"^ 

(31) -/ cos{y\A,\)e~"-^dy^-=T-'^^e-^. 
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This is a special of the Bessel integral identities known as Weber-Sonine formula. See for instance 
Andrews, Askey and Roy (1999, p. 222) and Watson (1966, p. 394 eq. (4)) for the identity and also 
those references for Bessel functions. It now follows rather immediately that, 

Proposition 4.1 For the model described by (6), let (ri, . . . , t„) have an arbitrary distribution where 
the joint Laplace transform has a known form. Then the marginal likelihood is given by, 



^(X|i?) 
where 



gnA/3 






-{y't/2+py2)T, 



Y[cos{yi\Ai\)dy^ 



E 



n 

-(y-'/2+/3V2)T. 



n 



e 
U=i 
is the joint Laplace transform o/ (ti, . . . , t„) evaluated at uji — yf /2 + /3^/2 for i — 1, . . . , n.O 

The next result which first appears in James (2005c) [see also James (2005b)], which can be 
thought of an unpublished earlier version of this manuscript, describes the case where Ti is repre- 
sentable as a functional of a Poisson random measure. Since positive Levy processes can be con- 
structed from Poisson random measures this represents a very rich class. 

Proposition 4.2 Let N denote a Poisson random measure on a Polish space S^ with sigma-finite 
mean intensity v, such that for each positive function g, the corresponding random variable N{g) 
has Levy exponent ^(wf?) = /^(l — e^^^^>'^)i'{dx). Suppose that Ti — N{gi) for positive functions 
(gi) on 3(j . Then since Y^^=i ^{^i9i) = ^(^7=i ^i9i)' *^ follows that (ti, . . . , r„) has the joint Levy 
exponent ^(fJ) — ^(X]r=i ^i9i)- Then for the model described by (6), the likelihood is given by, 

pnAp 



^(X|^) = -^ / e-''^''^]lcos{yM^\)dy^ 
where il(x) — X]"=i ^i9i{^) with uJi — yf /2 + 0^/2 for i — 1, . . . , n.D 
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Remark 16. Notice that we have stated the resuh in terms of quite arbitrary (ri, . . . , t„). This 
is because the expression (31) has nothing to do with the distributional properties of r. 

Remark 17. Hereafter we set 

n 

(32) '^{y\^i) = l[cos{y,\A,\) 

Remark 18. The appearance of integrals involving Besscl functions is certainly not new to 
applications in finance as can be seen in the case of the important work of Yor (1992) on Asian 
Options. See also Carr and Schroder (2004). 

5 General OU likelihoods 

We now apply Proposition 4.1, in the case of where (ri, . . . , r„) are based on the integrated OU mod- 
els described by (11). The task is to calculate the joint Laplace transform evaluated at (wi, . . . , u;„). 
This is straightforward from the construction given section 2.2. which implies that 

n n — 1 

A^W^Ti = Siw(O) + ^[sj+iO;n + \Zi - 0;]w;] + [Z„ - 0„]tJ„ 
i=\ l=\ 

where for / == 1, . . . , n, s; == (1 - cr^^)^l^i^i^^^'^^^^^% Then it is not difficult to see that the 
joint Laplace transform of (ri, . . . , t„) is of the form 



(33) Li(y|^)=e-'^^^^^e 



v{s^)^-^{^^) 



n-\ 

TT g-$(a;.|t..) 

«=1 



where terms are explicitly defined in the next result which gives the likelihood. 

Proposition 5.1 For the model described by (6), let (ri, . . . ,r„) be defined by the OU models as 
in (11). Then the marginal likelihood in (10) is, 



nAfJ r « 



where cOi — yf /2 + /3^/2 and Vi — r^Si+i. '^{y\^) is defined in (32) and Li(y|i?) is the joint Laplace 
transform evaluated at (a>i, . . . ,Cij„) with form specified by (33). The Levy exponents in (33) are 
specifically defined as follows, for a = XA, 

(i) ^LU,\v,) = /g_ ^b{X-^[v,u + u;,{l - u)])^, fori^l,...,n-l 

(m) <&(a;„) = jl_^ 7A(A-i^„(l - u))^ 

(iv) ip{si) — L ip{siX^^u) — , is the Levy exponent of v{Q) evaluated at siX^^. 
D 

Consider now the following result which we will return to in section 6. 
Proposition 5.2 Consider ^{uji\vi), $(wi), and let A{vi\LUi) = ^(uji\vi) — ^{oJi). Define 



Dp{y,e''y\u;,)= e-^^'p{ds)- e^"^^' p{ds) 

Jy JyC^ 



Then 
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(i) A(i;,|c^,) = a/o'/o"(l 



^'''')e-^^^^-''^'' p{ds)Fa{du). 



(ii) K{v,\LOi) = j^{l - e-^'y)Dp{y,e^y\w,)y-^e-ye^^ydy. 

(Hi) It follows that for each fixed uji, K{t\ijJi) is the Levy exponent, evaluated at t, of an infinitely 
divisible random variable with Levy density Dp{y, e"'y\wi)y^^ e^^ e^'^ . 

Proof. Note that A(wi|wi) = a jQ^^[ip{[viU + uJi{l-u)])-ip{[uJi{l-u)])]Fa{du). Statement (i) 
is simply the Levy density representation of this. Statement (ii) foUows by the change of variable 
y = MS, and exploiting the scale invariance of the measure u^^du. □ 



This allows one to better understand the representation of the joint Laplace transform 



(34) 



Li(y|^) ^c-'V{si)^-'S>(w^) 



TT ^-<I.(t^,|t,,) 



^-ipCsi) 



n- 



,-*(^0 



n-1 



-AKI^,) 



where Wj depends only on yi and each Vi depends on (y^+i, . . . , j/„). We also note that 



(35) 



L2(y|^):-c-'^(^i) 



.4=1 



~$(Wi 



is also a joint Laplace transform. In fact examining (35) more closely we see that it is the joint Laplace 
of a sequence random variables (^1, ... ,^„), where A^i — Civ{0) + [Z.i — Oi].lieTeci — (1— c^")c^"(*^^^ 
However note that if Ci — {I ~ e^"), then when v{t) is stationary, it follows that the marginal 
distribution of this version of ^i is equivalent to t^. Since we later propose the use of joint densities 
based on (34) and (35) one may want to vary the value of Ci in (35) as this may increase accuracy. 

6 Some distribution theory for OU-FGGC models 

We have already mentioned that the class of infinitely divisible random variables which are GGC's are 
closely linked with Dirichlet process mean functionals. When the Gamma process has a finite shape 
measure say 9H, then every such GGC can be expressed as TgMgn where Mgn = /n xPenidx) is 
a Dirichlet process. We will call such GGC's finite GGC's or FGGC. One implication is that one 
may apply some of the distribution theory we have developed for the OU-F to these models. In this 
section we shall assume that Z is derived from a finite GGC and demonstrate some nice properties 
of the corresponding OU process which are also relevant to sampling likelihoods and option pricing 
calculations. First note that if Z is an FGGC then its Levy density and Levy exponent are given by 



(36) 



V'(w) = / de{uJx)H{dx) and Oy^^ / c-y^''H{dr) 
Jo Jo 



where i7 is a probability measure 

Remark 19. The term finite GGC should not be confused with the term finite activity. That 
is to say FGGC are infinite activity models as can be seen from their Levy density in (36). 



6.1 Results for perfect sampling relevant OU-FGGC components 

Proposition 6.1 Suppose that Z is a BDLP with specifications given in (36). Then the following 
results hold. 
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(i) In the stationary case the corresponding OU process v{t) is such that v{0) is a non-finite GGC 
with Levy exponent 

dg{u!u)S{u)u~ du = / de{Luux)H{dx)u~ du 

Jo Jo 

where S{u) — J H{dy) is a survival function. 

(a) Consider the Levy exponent $(a>) described in Proposition 4-1- Then in this setting it takes 
the form 

(37) $H=/ dga{cjr)Qa{dr) = / dga{cj{l - u)x)H{dx)Fa{du) 

Jo Je-» Jo 

where Qa is a probability measure corresponding to the distribution of a random variable R — 
(1 — U)W where U has distribution Fa and W is independent of U and has distribution H. 

(Hi) Equivalently Z{a) — Yga = Jq (1 — e^^ y')Z{dXy) is a random variable with Levy expo- 
nent (37) and hence is a finite GGC and has the representation Z{a) — Yga — TgaMgaQ^ 

(iv) If the support of H is finite then the support of Qa is finite. This implies that the perfect 
simulation method described by (17) and (18) applies to MgaQ^- One may choose uMgaQ^^-N 
and lMeaQa.,-N , according to the upper and lower support points ofQa- Bn.ga is Beta {1,6a) 

and X„ = R = W (1 — U) has distribution Qa- 
(v) For 9a = 1 the density of Mq^ has the form, 

n 

where Qa{x) — L Qa{dt). Densities for 9a > 1 are obtained by substituting 9a — 9 and Qa = H 
in (16). 

(vi) The Levy exponent of Yga ~ L e^^^ ^i Z{d\y) is similar to (37) but with dgai^ux) in 
place of dga{uj{l — u)x). Hence Yga — TgaMg^Q^, where Qa corresponds to the distribution of 
R^WU. 

(vi) Yea converges in distribution to ve{0) as e^" -^ 0. 

D 



Proof. The proof of (i) and (ii) are obvious by substituting the form of ip in (36) into (8) 
and (9). The remaining results follow as consequences. The density in {v) is obtained from Cifarelli 
and Regazzini (1990) or Cifarelh and MehUi (2001). □ 

The next result is rather curious but as we shall show can play a powerful role in Monte Carlo 
procedures. 

Proposition 6.2 Consider the setting in Proposition 6.1 then the Levy exponent K{t\uji) described 
in Proposition 5.2 takes the form 

/•oo 

A(t|wi) = / dgaitr)Qai^^{dr). 
Jo 
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where Qa\uji is a probability measure corresponding to a random variable 

UW 



R 



1 + W{1 - U)uji 



where U has distribution Fa and W has distribution H. As a consequence results analogous to 
Proposition 6.1 apply to this setting. 

Proof. Similar to Proposition 5.2 wc examine 

K{vi\uji) ^ a I [tp{[viU + uji{l- u)]) -^'{[^ii'^- u)\)]Fa{du). 
However in this case ^p{[viU + LUi{l — u)]) — 'ip{[LUi{l — u)]) is equivalent to 



[dg{[v^u + u!i{l -u)]y) - de{[Lu^{l -u)]y)]H{dy) 

'0 

Now using properties of the natural logarithm it follows that 

dedviu + Wj(l - u)]y) -~ ^^([^^(l - u)]y) = do — '- — 

VI + tjj(l -u)yj 

concluding the result. D 

6.2 OU-FGGC Monte Carlo Densities 

The next result, whose present importance is that it can be used effectively in Monte Carlo sim- 
ulation procedures, follows immediately from Propositions 6.1 and 6.2 and standard augmentation 
arguments. 

Theorem 6.1 Suppose that the joint Laplace transforms Li(y|'!9), and hence L2(y|'!9), satisfies the 
conditions in Proposition 6.1 and 6.2. Specify uji ~ {yf + /3^)/2. From this, for i ~ 1, . . . ,n, we can 
let (Ti,Mi) = {Tga,TeaMgaQa) denote iid pairs of random variables. Similarly, independent of the 
above sequence, define independent pairs (Gi, M^.), where Gi = Tga o-nd independent of Gi, M^^ 
has the distribution of a mean functional described in Proposition 6.2 for fixed uJi. Let 

Si = (w(0),(T„M,),(G„M^J) 

denote the joint vector ofAn + 1 independent components, with joint density fsi{- I'd, y). Similarly let 
■=■2 = (v(O)j {Ti, Mi)) denote the joint vector of 2n+l independent components with density /h^ ('I'?) 
specified by Proposition 6.1 and not depending on y. Then, 

(i) Li(y|^) -E[e--i-(0)]nr=iE[e-'^'^*^inr=iE[e""'''''^"'] 
(ii) L2(y|i?) =E[e--i-(o)]nr=iE[e-"'^*^'] 

(Hi) Suppose that L„ Li(y|^)ni=i '^Vi ^ °°' then by augmenting the expression in (i) there exists 
a joint density o/(Si,Y) given by 



n n 



(38) /H,(Ci,y|^)« e-^^''ne-"'*'"-n^"''"'""'/s,(Cil^,y) 

where Ci = (v, (ti,mi), {ui,m^.)), with obvious meaning. 
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(iv) Suppose that Jjg„ L2(y|'*^) 0"=! '^Vi '^ °°' then by augmenting the expression in (ii) there exists 
a joint density of (S2 , Y) given by 



/H.(C2,y|^) ex e-^^'^n^""'*'"''/-^(C2|^,y), 



where C2 = (w, [U^mi)). 



(v) Writing w(0) — TgMg and integrating out all the Gamma random variables in (Hi) it follows 
that there exist a joint density of {Mg, (Mi), {M^.),'Y) given proportional to 

n 

(39) /M,(i)(l + sity' n (1 + ^^m^^'"{^ + v^r,)-'"'fMArn^)fM^^ in). 

i=l 

6.3 OU-FGGC option pricing densities 

The last result, extends Proposition 3.4 and again is pertinent to the option pricing formula discussed 
in BNS(2001a, 6.2) and Nicolato and Vernados (2003). 

Proposition 6.3 Let x*{t) be defined by the BDLP Z which is an FGGC with specifications (36). 
Additionally, for < s < t, set A = {t — s) and define h{A,s) = (1 — e^^ )v{s) and /i* = 
liA + x*{s) +/3/i(A,s). Then the conditional density of x* (t)\x*{s),v{s) is given by 

(f>{x\nl + Py, h{A, s) + y)q0a{y)dy 

where qea{y) ~ L '^9a{y\v)fMgaQ^{^)d''^- With the density further described by the specifications in 
Proposition 6. in 

7 Some practical issues for general OU likelihood estimation 

The likelihoods given in Propositions 4.1, 4.2 and 5.1 serve the purpose of integrating out the infinite- 
dimensional nuisance parameters. A natural question is how to exploit these results in a practical 
sense. In the forthcoming sections we shall focus on the OU models but many parts of our discussion 
can be extended to more general processes where the joint Laplace transform has an accessible form. 
Our goal at minimum will be to discuss ways to evaluate the likelihood by Monte Carlo procedures. 
This could then be used in conjunction with simulated maximum likelihood estimation or other 
such techniques. Similar to section 3.4 we will also be thinking about Bayesian type estimation 
procedures. That is, we wish to calculate 



(40) E[/i(t?)|X] = 



/^ h{d)n{dd)e-^pu{y\^)^{y\ii) nr=i % 



/^ ^{d^)e^APL^ {yW^{y\^l) WU dy^ 
where we set =^ = (M!^, 0). 

7.1 Calculating Levy exponents 

In order to utilize Proposition 5.1 one needs a manageable expression for 

/■I ri r^ 

(41) $(u;i|w2)=/ ^[L02U + L0i{l-u)])—-~a I ^{[lo2U + uji{l - u)])Fa{du) 

Je^a U jQ-a 
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where loi and lu2 just denote two arbitrary non-negative numbers. We have removed the dependence 
on the scale factor A^^, which can otherwise be absorbed in {101,102)- We wiU assume that ip has a 
known form. One can see that (41) is the Levy exponent of the joint distribution of (7). However, 
even if we wished to try to apply a direct inversion the results for the OU-F would suggest that, in 
general, the joint density of (7) has a rather non-obvious form. That is to say, except for the OU-F 
case, it is probably just as well to work directly with (41). Now again note importantly that in order 
to calculate (41) we only need knowledge of tp and not the Levy density p of Z. In many cases 
manual evaluation of (41) may not be obvious. One can then resort to numerical methods available 
in standard mathematical packages or one can carry out a one time Monte-Carlo approximation 
based on the following, somewhat obvious, result. 



Proposition 7.1 Let Ui for I = 

B^L0i\t02) =Ef=lH'^2Ul + iOi{l 



D 



1,...,B denote iid random variables with distribution Fa 
-Ui)) Then 



Let 



E 



$(u;i|tJ2) =$(wi|w2) 



Note that our intention is to use a one time calculation of ^{loi\lo2), based on large i?, to get a highly 
accurate approximation to <^{loi\lo2)- Our intent is not to continuously generate different realizations 
of $(a;i|aj2) within a loop. In other words one stores a set of (C/;). The remaining sections will assume 
that we have been able to get an expression for <^{loi\lo2) by some means. 

Remark 20. As seen from our results in section 6 we do not necessarily need to work with (41) 
in the case of OU-FGGC models. 



7.2 Monte Carlo method 

It is well-known that classical iid Monte-Carlo, MCMC and SIS procedures are well-suited to high- 
dimensional integrals. However, at first glance, one might think it is difficult to work with the 
expressions involving cosines. Specifically our likelihoods are expressed in terms of '^{y\p) which 
oscillates between positive and negative values. On the other hand, we note that 

|<^(y|M)|<|cos(yi|Ai|)|<l 

for all (yi, . . . , j/„), which suggests that a product of cosines is not any more unstable than a single 
cosine. Monte Carlo procedures just require a reasonable proposal density and otherwise deal with 
terms such "^(yj/x) in terms of an expectation E[ft,(Y)] where h depends on ^(y|^) and possibly other 
terms. Accuracy then becomes primarily a function of the number B of computer iterations. That 
is, in terms of B Monte Carlo replications. This is in contrast to numerical techniques which have 
difficulty handling high dimensions in n. See for instance Liu (2001), Chen, Shao and Ibrahim (2000) 
and Kong, Liu and Wong (1997). 

The idea of Monte Carlo in the general setting is in principle no different than that outlined in 
section 3.4. Except now we will sample from densities built from Li(y|'i?) and L2(y|'i?). Similar to 
Theorem 6.1, this would be possible if its prospective normalizing constant was finite. Note one can 
sample these densities without explicit knowledge of the normalizing constant via MCMC methods. 
We now give a description of its normalizing constant. 

Proposition 7.2 Let £_i — Civ{0) + [Zi — Oi] , for i — 1, . . . ,n. Then max(cii;(0), [Z^ — O^]) < ^i < t^, 
and the following results hold. 



(i) N^,i = u Li(yi^) nr=i dy^ = ^"E nr=i 



T" e 



-/3" 



(n) N^.2 = /k„ L2(y|^) nr=i dy^ :- ^"E [iTU S? 



2TTTi 



<7r"E 
< 7r"E 



n" ^- 
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Proposition 7.2 follows by a straightforward argument which can be seen more clearly in section 8. 
We see from Proposition 7.2 that the prospective normalizing constants are just based on negative 
moments of the random variables r^, ^i, w(0) and Zi — O^ which may or may not exist. One can 
always ensure finiteness by adding a small positive constant to any of the random variables. For 
instance one uses the model based on (ti + 6, . . . , t„ + b) for a small 6 > 0. Hereafter we shall then 
assume that modification is made if deemed necessary. This now allows us to describe two possible 
densities for Monte Carlo implementation as follows 

(42) N^,igi(y|^) = Li(y|^) 
and 

(43) N^,2Q2(y|^) = L2(y|^). 

Naturally, from the point of view of Monte Carlo (theoretical) accuracy, Qi{y\^) is the most 
desirable. However, Q2{y\'d) is in general easier to sample from. One can also adjust Q2{y\'d) further 
if necessary. Define the ratio 

Proposition 7.3 Consider the densities defined in (4-2) and (4-3) and the Bayesian posterior quan- 
tity given in (40). Additionally lefE^j denote expectation with respect to the respective joint density 
ofY = (Yi, . . . ,Yn), Qj{y\'&) for j = 1, 2. Define also Ej to denote expectation with respect to the 
joint densities 7r(i9)(5j(y|i?) forj = 1,2 Then it follows that 

(i) J^{^\^)^Ko^^^^.l['^{Y\^i)] 

(ii) ^(X|i?) =iV^,i^E^,2[^(Y|M)*(Y|^)] 
(Hi) This implies that 

^ ' ^ ^ >\ ^ Ei[e"^/37V^,i'^(Y|M)] E2[e"^'37V^,2'^(Y|M)*(Y|^)] 

Hence a Bayesian approach proceeds similar to section 3.4 by sampling (li,;, . . . ,y„,;,i?/) for I = 
1, . . . B times from either '^{'d)Qi{y\'d) or '^{'&)Qi{y\'&) and put them into appropriate empirical 
versions of (44). We now say a few more words about sampling from the respective densities 

7.2.1 Sampling from Qi 

In general an exact expression for the conditional marginals of say Yk\Yi, . . . ,Yn based on Qi{y\'&) 
can be worked out but it is a bit tricky. As such we do not discuss this. Note that for OU-FGGC 
models one one can definitely use Theorem 6.1 to sample from the joint distribution of (Si, Y,?9) 
based on (38) or sampling based on the density (39). These methods are facilitated by the fact that 
we can use the perfect simulation methods described in section 2.3.2., with specifications given by 
Proposition 6.1 and 6.2. 

7.2.2 Sampling from Q2 

Notice that in general Q2(y|^) has an almost independent structure and hence a rejection sampling 
procedure is straightforward. If however we know the distribution of vq we can introduce a further 
augmentation based on 

/■oo 
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where again si = Y^"^i Ci{yf + l3'^)/2. Hence the Monte Carlo procedure can be based on a joint 
density of (Y, y|^)|^ given as 



j(y,w) « 



n« 



-y,vCi/2 -i(u,i) 



-/3=Er..c./2. („) 



In the OU-FGGC case we may again use Theorem 6.1 in an obvious way. 

8 General approach 

So far we have advocated the idea of samphng using the joint Laplace transform or some variation 
of that. Since we focused on the BNS models we were able to highlight some nice features. However 
our claim is that one can implement similar procedures. This leads us to derive a similar approach 
that is influenced by some arguments in Devroye (1986a) but where we do not necessarily sample 
using the Laplace transform. That is we give another representation of the likelihood that can be 
numerically evaluated via the simulation of random variables. First let p = (pi, . . . ,p„) denote a 
vector of positive numbers and for each i, let 



H{yi\pt) 



zC ^"i for j/i > 



\/2np. 
denote a half Normal density. Now, notice that < 1 — Y[i=i cos(2/i) < 2, and 



(45) 



l-Y[cos{y,\A,\) 



H{yi\pi)dyi = 1 -e 



Cn{A,p) 



This follows from applications of the Fourier-Cosine identity that we used in section 4.1. From these 
facts we describe a joint density 

Proposition 8.1 Augmenting the expression in (4-5) leads to a joint density of an array of positive 
random variables Y — {Yi,n, ■ ■ ■ , Yn.n} given by, 



rn{y\p) 



[I -nUcos{y,\A\)]K=iH{yM) 



C„(A,p) 

Equivalently, for k — l,...,n, the conditional density o/ Yfe.„|Yi. „,..., Yfe-i^n is proportional to 
[1- XkCos{yk\Ak\)]H{yk\pk), where Xk^ e ^^=k+i 2 JJ^^.^ cos{yi\A,\) for k = 2, . . . ,n - 1, \i = 



e ^?=2 2 \ and A„ = Y["=i cos(yj|Ai|). 

Define the function, verified via Fubini's theorem and standard Normal integration. 



T„(^) := i 



E 



n- 



-(y.V2+/3V2)T. 



Yldy^^E 



n 



.-/3V, 



^ V27rTi 



< 



n:i 



These points lead to following representation of the likelihood. 



Proposition 8.2 Suppose that for fixed n, 
4-1 may be written as 



llj=l ^/T- 



< 00, then the likelihood in Proposition 



^Afi 



T„(^)-^^^^%^E[f7(yi^„,...,y„, 
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where 



E 
f^(2/i,...,y„|^)-- 



n: 



-fe.V2+r/2)n 



m=lH{y^\p^) 

and the random vector {Yi,„, . . . , Yn_n} has its joint distribution described by Proposition 8.1.0 

Remark 21. Proposition 8.2 shows that one may approximate the hkchhood by simulating 
random variables described in Proposition 8.1. Such an approach should work well with a Bayesian 
procedure. Methods to easily sample the random variables in Proposition 8.1, may be deduced from 
Devroye (1986a, b). In fact, through a personal communication with Luc Devroye we were informed 
that one at time sampling using the conditional distributions in Proposition 8.1 is routine as it 
constitutes essentially a sampling from a Normal density times a factor between and 2. Hence 
rejection sampling is easy and furthermore the normalizing factor is not needed. One may also use 
other densities. 

Remark 22. Note that one needs also to evaluate T„(?9). Of course this can also be done by a 
Monte Carlo procedure using the density in Proposition 8.1. 

9 Examples 

In this section we will present some examples where we sketch out a few details related to our 
exposition. We will not concern ourselves too much with constants. Note that all the examples 
presented are infinite- activity processes. In the case where the distribution of w(0) is not obvious we 
would simply approximate it when it is based on OU-FGGC models using Proposition 6.1, or choose 
an arbitrary law for w(0) in a more general setting. 

9.1 OU-Stable 

Suppose that Z is stable subordinator of index a specified by ip{i^) = tj". Then it is known, or oth- 
erwise obvious, that v{t) also has a stable law of index < a < 1 with Levy exponent a>" /„ u°'~^du. 
Notice that the Levy exponent of the corresponding 

$(w) = w'^ /" (1 - ufu-^du 
Jc-° 

corresponds also to a stable law of index a. Here, for simplicity of presentation, suppressing constants 
and setting [3 — Q we may use Q2 which is based on sampling the joint Laplace transform 

n 

(46) e-E"=i^''l°[|e-^''° 

Noting the simplicity of (46) it is good to recall that in general the densities of a stable law are 
only known in a complicated form. So here is a case where a Laplace transform approach is perhaps 
preferable despite the availability of the relevant densities. A nice exception to the preceding comment 
is when a = 1/2 corresponding to an inverse Gamma law of index a ~ 1/2. However in that case 
(46) is 

n 

For further simplification wc may use the augmentation procedure described in section 7.2.2 applied 
to (46) to get 



J|e-^"e-"^'/„(«) 



i=l 
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where /„ corresponds to a stable density. Note that although the stable density can be complicated 
there are many routines available to easily sample stable random variables. 

Remark 23. The Stable law process produces a log price process with heavy tails which may 
not be desirable for all applications. However see the work of Carr and Wu (2003). Additionally, we 
note that it would not be tremendously difficult to use Qi in this case. 

9.2 IG-OU 

This example is based on the calculations given in Barndorff-Nielsen and Shephard (2003, p. 292) 
where v{t) has an Inverse Gaussian distribution. Here, letting Ci,C2 denote constants and setting 
/3 = 0, by BNS(2003, eq. (54)) one has 

$(u;) = -y^C\ f (1 - u)-\i{l + Chy^uy'^^du. 



BNS(2003) show that this can be written in terms of the hyperbolic arc-tangent function[see also 
Nicolato and Vernardos (2001) and Carr, Geman, Madan and Yor (2003)], we do not repeat that 
here. Note however by using the fact that v{t) has a Inverse Gaussian distribution one can work 
with the augmented version of Q2 which is proportional to 



3/2p-K7't'+<5'«"'] TTp-y^p-*K) 






for appropriate values of 7 and 5 and is not difficult to sample from. 

9.3 OU-LogNormal 

Suppose that Z is based on a LogNormal distribution with density 

f[x) = ^ic-^(i°s(^»' for X > 0. 
\/2tt X 



We have chosen this example because, despite the fact that it has a density with a nice closed form, 
its corresponding Levy density p is unknown. Despite this we can still use a sampler based on Q2. 
This is because its Levy exponent is given by 



■(/'(w) = — log 



00 1 1 



2tt X 



This can be numerically approximated hence the relevant quantities $(a;) can then be numerically 
approximated. Again this approximation should be done before the main Monte-Carlo procedure is 
used. Note that one would find it difficult or impossible to employ a series approximation in this 
case, as it depends on knowledge of p. 

9.4 OU-FGGC where H is the Arcsine distribution 

Here we close with one of the more interesting examples of known FGGC models. In this setting let 
H be the Arcsine law, that is there is a corresponding random variable W which is Beta(1/2, 1/2). 
Cifarelh and Melilh (2000) show that in this setting for all 6* > 0, Men is Beta(6' + 1/2,61 + 1/2). 

Hence Z[t) = TgtB(^gt+i/2.et+i/2), where here -B(6/t+i/2,ei+i/2) means a beta random variable with 
parameters indicated in the subscript. In this case the distribution oi R = (l — U)W described in 
Proposition 6.1 has bounded support on [0, 1]. Hence we may apply the perfect sampler both for 
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option pricing and Monte Carlo methods specified according to Theorem 6.1 and Propositions 6.1 
and 6.2. That is apply section 2.3.2 to sample from (Si, Y). To be clear given Yi one may draw M^^. 
from /mu,- by using Proposition 6.2 and creating uMga.-N = 1 a-nd lMga.wi,-N = 0, B^fia is Beta 

(1, Oa) and Xn ~ UW/{1 + W{1 — U)wi). Then one draws W from the Arscine law and U from Fa 
to get X(n). Draws from more complex densities for Mu,. can then be obtained by other standard 
methods. One can also work with the exact form of the densities via Cifarelli and Regazzini (1990). 
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